Load libraries.

library(Seurat)
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio
Attaching SeuratObject
library(URD)
Loading required package: ggplot2
Loading required package: Matrix
Registered S3 method overwritten by 'gplots':
  method         from 
  reorder.factor gdata
library(tidyverse)
Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.2 ──✔ tibble  3.1.8      ✔ dplyr   1.0.10
✔ tidyr   1.2.1      ✔ stringr 1.5.0 
✔ readr   2.1.3      ✔ forcats 0.5.2 
✔ purrr   1.0.1      ── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ tidyr::expand() masks Matrix::expand()
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
✖ tidyr::pack()   masks Matrix::pack()
✖ tidyr::unpack() masks Matrix::unpack()

Load merged eb, lb, ehf.

load("/Volumes/Bud_SSD/20230122-urd/merged-seurat.Robj")

Grab count data and meta data.

count.merged <- as.matrix(merged.seurat@assays$RNA@counts)
Warning: sparse->dense coercion: allocating vector of size 6.2 GiB
meta.merged <- as.data.frame(merged.seurat@meta.data)

rm(merged.seurat)

Create urd.

merged.urd <- createURD(count.data = count.merged, meta = meta.merged, min.cells = 3, min.counts = 0)
2023-01-23 07:20:42: Filtering cells by number of genes.
2023-01-23 07:20:59: Filtering genes by number of cells.
2023-01-23 07:21:23: Filtering genes by number of counts across entire data.
2023-01-23 07:21:46: Filtering genes by maximum observed expression.
2023-01-23 07:22:11: Creating URD object.
2023-01-23 07:22:15: Determining normalization factors.
Warning: sparse->dense coercion: allocating vector of size 6.0 GiB2023-01-23 07:22:28: Normalizing and log-transforming the data.
as(<dgeMatrix>, "dgCMatrix") is deprecated since Matrix 1.5-0; do as(., "CsparseMatrix") instead
2023-01-23 07:23:12: Finishing setup of the URD object.
Warning: sparse->dense coercion: allocating vector of size 6.0 GiB2023-01-23 07:23:26: All done.

Add stage.

merged.urd@group.ids$stage <- as.character(merged.urd@meta[rownames(merged.urd@group.ids),"stage"])

Add embryo.

merged.urd@group.ids$embryo <- as.character(merged.urd@meta[rownames(merged.urd@group.ids),"embryo"])
#embryo <- sort(unique(merged.urd@group.ids$embryo))
embryo <- c("EB1", "EB2", "AVJ4-56", "AYQ8-11", "AYR4-57", "BBH6-17", "AXL7-21", "EHF1", "EHF2")

Calculate variable genes. Use only Flox embryos.

var.by.embryo <- lapply(c(1, 2, 3, 4, 5, 6, 7, 8, 9), function(n) {
  findVariableGenes(merged.urd, cells.fit=cellsInCluster(merged.urd, "embryo", embryo[n]), set.object.var.genes=F, diffCV.cutoff=0.3, mean.min=.005, mean.max=100, main.use=paste0("Embryo ", embryo[n]), do.plot=T)
})

Add variable genes.

var.genes <- sort(unique(unlist(var.by.embryo)))
merged.urd@var.genes <- var.genes

Calculate PCA.

merged.urd <- calcPCA(merged.urd, mp.factor = 2)
[1] "2023-01-23 07:32:26: Centering and scaling data."
[1] "2023-01-23 07:32:34: Removing genes with no variation."
[1] "2023-01-23 07:32:36: Calculating PCA."
[1] "2023-01-23 07:37:28: Estimating significant PCs."
[1] "Marchenko-Pastur eigenvalue null upper bound: 1.40724420496983"
[1] "31 PCs have eigenvalues larger than 2 times null upper bound."
[1] "Storing 62 PCs."
pcSDPlot(merged.urd)


set.seed(19)
merged.urd <- calcTsne(object = merged.urd)
plotDim(merged.urd, "stage", plot.title = "tSNE: Stage")
Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as of ggplot2 3.3.4.

plotDim(merged.urd, "cluster.whole", plot.title = "tSNE: Cluster")

Save.

save(merged.urd, file = "merged-tsne-urd.Robj")

Load.

load("merged-tsne-urd.Robj")

Calculate diffusion map.

merged.urd <- calcDM(merged.urd, knn = sqrt(dim(merged.urd@meta)[1]), sigma=NULL)
[1] "destiny determined an optimal global sigma of 41.63"
Warning: You have 1603 genes. Consider passing e.g. n_pcs = 50 to speed up computation.as(<dsCMatrix>, "dgTMatrix") is deprecated since Matrix 1.5-0; do as(as(., "generalMatrix"), "TsparseMatrix") instead
plotDimArray(merged.urd, reduction.use = "dm", dims.to.plot = 1:8, outer.title = "Diffusion Map (Sigma NULL, 214 NNs): Stage", label="stage", plot.title="", legend=F)
Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as of ggplot2 3.3.4.

plotDim(merged.urd, "stage", transitions.plot = 10000, plot.title="Developmental stage (with transitions)")


plotDim(merged.urd, "cluster.whole", transitions.plot = 10000, plot.title="Developmental stage (with transitions)")

Add cluster and stage.

merged.urd@meta <- merged.urd@meta %>% unite("cluster.stage", cluster.whole:stage, remove = FALSE)
merged.urd@group.ids$cluster.stage <- as.character(merged.urd@meta[rownames(merged.urd@group.ids),"cluster.stage"])

Calculate pseudotime.


pseudotimePlotStabilityOverall(merged.urd)

plotDim(merged.urd, "pseudotime")


plotDists(merged.urd, "pseudotime", "stage", plot.title="Pseudotime by stage")


plotDists(merged.urd, "pseudotime", "embryo", plot.title="Pseudotime by embryo")


plotDists(merged.urd, "pseudotime", "genotype", plot.title="Pseudotime by genotype")

Save.

save(merged.urd, file = "merged-pseudotime-urd.Robj")

Extract pseudotime.

pseudo.data <- data.frame(rownames(merged.urd@pseudotime), merged.urd@meta$stage, merged.urd@meta$embryo, merged.urd@meta$genotype, merged.urd@meta$cluster.whole, merged.urd@meta$germ.whole, merged.urd@pseudotime$pseudotime)

save(pseudo.data, file = "pseudotime-urd.Robj")
LS0tCnRpdGxlOiAiVVJEIgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCgpMb2FkIGxpYnJhcmllcy4KYGBge3J9CmxpYnJhcnkoU2V1cmF0KQpsaWJyYXJ5KFVSRCkKbGlicmFyeSh0aWR5dmVyc2UpCmBgYAoKTG9hZCBtZXJnZWQgZWIsIGxiLCBlaGYuCmBgYHtyfQpsb2FkKCIvVm9sdW1lcy9CdWRfU1NELzIwMjMwMTIyLXVyZC9tZXJnZWQtc2V1cmF0LlJvYmoiKQpgYGAKCkdyYWIgY291bnQgZGF0YSBhbmQgbWV0YSBkYXRhLgpgYGB7cn0KY291bnQubWVyZ2VkIDwtIGFzLm1hdHJpeChtZXJnZWQuc2V1cmF0QGFzc2F5cyRSTkFAY291bnRzKQptZXRhLm1lcmdlZCA8LSBhcy5kYXRhLmZyYW1lKG1lcmdlZC5zZXVyYXRAbWV0YS5kYXRhKQoKcm0obWVyZ2VkLnNldXJhdCkKYGBgCgpDcmVhdGUgdXJkLgpgYGB7cn0KbWVyZ2VkLnVyZCA8LSBjcmVhdGVVUkQoY291bnQuZGF0YSA9IGNvdW50Lm1lcmdlZCwgbWV0YSA9IG1ldGEubWVyZ2VkLCBtaW4uY2VsbHMgPSAzLCBtaW4uY291bnRzID0gMCkKYGBgCgpBZGQgc3RhZ2UuCmBgYHtyfQptZXJnZWQudXJkQGdyb3VwLmlkcyRzdGFnZSA8LSBhcy5jaGFyYWN0ZXIobWVyZ2VkLnVyZEBtZXRhW3Jvd25hbWVzKG1lcmdlZC51cmRAZ3JvdXAuaWRzKSwic3RhZ2UiXSkKYGBgCgpBZGQgZW1icnlvLgpgYGB7cn0KbWVyZ2VkLnVyZEBncm91cC5pZHMkZW1icnlvIDwtIGFzLmNoYXJhY3RlcihtZXJnZWQudXJkQG1ldGFbcm93bmFtZXMobWVyZ2VkLnVyZEBncm91cC5pZHMpLCJlbWJyeW8iXSkKYGBgCgoKYGBge3J9CiNlbWJyeW8gPC0gc29ydCh1bmlxdWUobWVyZ2VkLnVyZEBncm91cC5pZHMkZW1icnlvKSkKZW1icnlvIDwtIGMoIkVCMSIsICJFQjIiLCAiQVZKNC01NiIsICJBWVE4LTExIiwgIkFZUjQtNTciLCAiQkJINi0xNyIsICJBWEw3LTIxIiwgIkVIRjEiLCAiRUhGMiIpCmBgYAoKQ2FsY3VsYXRlIHZhcmlhYmxlIGdlbmVzLiBVc2Ugb25seSBGbG94IGVtYnJ5b3MuCmBgYHtyfQp2YXIuYnkuZW1icnlvIDwtIGxhcHBseShjKDEsIDIsIDMsIDQsIDUsIDYsIDcsIDgsIDkpLCBmdW5jdGlvbihuKSB7CiAgZmluZFZhcmlhYmxlR2VuZXMobWVyZ2VkLnVyZCwgY2VsbHMuZml0PWNlbGxzSW5DbHVzdGVyKG1lcmdlZC51cmQsICJlbWJyeW8iLCBlbWJyeW9bbl0pLCBzZXQub2JqZWN0LnZhci5nZW5lcz1GLCBkaWZmQ1YuY3V0b2ZmPTAuMywgbWVhbi5taW49LjAwNSwgbWVhbi5tYXg9MTAwLCBtYWluLnVzZT1wYXN0ZTAoIkVtYnJ5byAiLCBlbWJyeW9bbl0pLCBkby5wbG90PVQpCn0pCmBgYAoKQWRkIHZhcmlhYmxlIGdlbmVzLgpgYGB7cn0KdmFyLmdlbmVzIDwtIHNvcnQodW5pcXVlKHVubGlzdCh2YXIuYnkuZW1icnlvKSkpCm1lcmdlZC51cmRAdmFyLmdlbmVzIDwtIHZhci5nZW5lcwpgYGAKCkNhbGN1bGF0ZSBQQ0EuCmBgYHtyfQptZXJnZWQudXJkIDwtIGNhbGNQQ0EobWVyZ2VkLnVyZCwgbXAuZmFjdG9yID0gMikKcGNTRFBsb3QobWVyZ2VkLnVyZCkKCnNldC5zZWVkKDE5KQptZXJnZWQudXJkIDwtIGNhbGNUc25lKG9iamVjdCA9IG1lcmdlZC51cmQpCnBsb3REaW0obWVyZ2VkLnVyZCwgInN0YWdlIiwgcGxvdC50aXRsZSA9ICJ0U05FOiBTdGFnZSIpCnBsb3REaW0obWVyZ2VkLnVyZCwgImNsdXN0ZXIud2hvbGUiLCBwbG90LnRpdGxlID0gInRTTkU6IENsdXN0ZXIiKQpgYGAKClNhdmUuCmBgYHtyfQpzYXZlKG1lcmdlZC51cmQsIGZpbGUgPSAibWVyZ2VkLXRzbmUtdXJkLlJvYmoiKQpgYGAKCkxvYWQuCmBgYHtyfQpsb2FkKCJtZXJnZWQtdHNuZS11cmQuUm9iaiIpCmBgYAoKQ2FsY3VsYXRlIGRpZmZ1c2lvbiBtYXAuCmBgYHtyfQptZXJnZWQudXJkIDwtIGNhbGNETShtZXJnZWQudXJkLCBrbm4gPSBzcXJ0KGRpbShtZXJnZWQudXJkQG1ldGEpWzFdKSwgc2lnbWE9TlVMTCkKCnBsb3REaW1BcnJheShtZXJnZWQudXJkLCByZWR1Y3Rpb24udXNlID0gImRtIiwgZGltcy50by5wbG90ID0gMTo4LCBvdXRlci50aXRsZSA9ICJEaWZmdXNpb24gTWFwIChTaWdtYSBOVUxMLCAyMTQgTk5zKTogU3RhZ2UiLCBsYWJlbD0ic3RhZ2UiLCBwbG90LnRpdGxlPSIiLCBsZWdlbmQ9RikKCnBsb3REaW0obWVyZ2VkLnVyZCwgInN0YWdlIiwgdHJhbnNpdGlvbnMucGxvdCA9IDEwMDAwLCBwbG90LnRpdGxlPSJEZXZlbG9wbWVudGFsIHN0YWdlICh3aXRoIHRyYW5zaXRpb25zKSIpCgpwbG90RGltKG1lcmdlZC51cmQsICJjbHVzdGVyLndob2xlIiwgdHJhbnNpdGlvbnMucGxvdCA9IDEwMDAwLCBwbG90LnRpdGxlPSJEZXZlbG9wbWVudGFsIHN0YWdlICh3aXRoIHRyYW5zaXRpb25zKSIpCmBgYAoKQWRkIGNsdXN0ZXIgYW5kIHN0YWdlLgpgYGB7cn0KbWVyZ2VkLnVyZEBtZXRhIDwtIG1lcmdlZC51cmRAbWV0YSAlPiUgdW5pdGUoImNsdXN0ZXIuc3RhZ2UiLCBjbHVzdGVyLndob2xlOnN0YWdlLCByZW1vdmUgPSBGQUxTRSkKbWVyZ2VkLnVyZEBncm91cC5pZHMkY2x1c3Rlci5zdGFnZSA8LSBhcy5jaGFyYWN0ZXIobWVyZ2VkLnVyZEBtZXRhW3Jvd25hbWVzKG1lcmdlZC51cmRAZ3JvdXAuaWRzKSwiY2x1c3Rlci5zdGFnZSJdKQpgYGAKCkNhbGN1bGF0ZSBwc2V1ZG90aW1lLgpgYGB7cn0Kcm9vdC5jZWxscyA8LSBjZWxsc0luQ2x1c3RlcihtZXJnZWQudXJkLCAiY2x1c3Rlci5zdGFnZSIsICIxLjFfZWIiKQoKbWVyZ2VkLmZsb29kcyA8LSBmbG9vZFBzZXVkb3RpbWUobWVyZ2VkLnVyZCwgcm9vdC5jZWxscyA9IHJvb3QuY2VsbHMsIG49NTAsIG1pbmltdW0uY2VsbHMuZmxvb2RlZCA9IDIsIHZlcmJvc2U9RikKCm1lcmdlZC51cmQgPC0gZmxvb2RQc2V1ZG90aW1lUHJvY2VzcyhtZXJnZWQudXJkLCBtZXJnZWQuZmxvb2RzLCBmbG9vZHMubmFtZT0icHNldWRvdGltZSIpCgpwc2V1ZG90aW1lUGxvdFN0YWJpbGl0eU92ZXJhbGwobWVyZ2VkLnVyZCkKCnBsb3REaW0obWVyZ2VkLnVyZCwgInBzZXVkb3RpbWUiKQoKcGxvdERpc3RzKG1lcmdlZC51cmQsICJwc2V1ZG90aW1lIiwgInN0YWdlIiwgcGxvdC50aXRsZT0iUHNldWRvdGltZSBieSBzdGFnZSIpCgpwbG90RGlzdHMobWVyZ2VkLnVyZCwgInBzZXVkb3RpbWUiLCAiZW1icnlvIiwgcGxvdC50aXRsZT0iUHNldWRvdGltZSBieSBlbWJyeW8iKQoKcGxvdERpc3RzKG1lcmdlZC51cmQsICJwc2V1ZG90aW1lIiwgImdlbm90eXBlIiwgcGxvdC50aXRsZT0iUHNldWRvdGltZSBieSBnZW5vdHlwZSIpCmBgYApTYXZlLgpgYGB7cn0Kc2F2ZShtZXJnZWQudXJkLCBmaWxlID0gIm1lcmdlZC1wc2V1ZG90aW1lLXVyZC5Sb2JqIikKYGBgCgpFeHRyYWN0IHBzZXVkb3RpbWUuCmBgYHtyfQpwc2V1ZG8uZGF0YSA8LSBkYXRhLmZyYW1lKHJvd25hbWVzKG1lcmdlZC51cmRAcHNldWRvdGltZSksIG1lcmdlZC51cmRAbWV0YSRzdGFnZSwgbWVyZ2VkLnVyZEBtZXRhJGVtYnJ5bywgbWVyZ2VkLnVyZEBtZXRhJGdlbm90eXBlLCBtZXJnZWQudXJkQG1ldGEkY2x1c3Rlci53aG9sZSwgbWVyZ2VkLnVyZEBtZXRhJGdlcm0ud2hvbGUsIG1lcmdlZC51cmRAcHNldWRvdGltZSRwc2V1ZG90aW1lKQoKc2F2ZShwc2V1ZG8uZGF0YSwgZmlsZSA9ICJwc2V1ZG90aW1lLXVyZC5Sb2JqIikKd3JpdGUuY3N2KHBzZXVkby5kYXRhLCBmaWxlID0gInBzZXVkb3RpbWUtdXJkLmNzdiIpCmBgYAoK